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Abstract 

We  extend  the  development  of  collocation  methods  within  the  framework  of  Iso¬ 
geometric  Analysis  (IGA)  to  multi-patch  NURBS  configurations,  various  boundary 
and  patch  interface  conditions,  and  explicit  dynamic  analysis.  The  methods  devel¬ 
oped  are  higher-order  accurate,  stable  with  no  hourglass  modes,  and  efficient  in 
that  they  require  a  minimum  number  of  quadrature  evaluations.  The  combination 
of  these  attributes  has  not  been  obtained  previously  within  standard  finite  element 
analysis. 

Key  words:  Isogeometric  analysis;  collocation  methods;  B-splines;  NURBS; 
explicit  dynamics. 


1  Introduction 


There  are  many  application  areas  of  finite  element  analysis  in  which  the  ef¬ 
ficiency  and  success  of  the  methodology  is  directly  related  to  the  number  of 
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quadrature  points  needed  to  integrate  arrays.  The  most  salient  example  is 
explicit  dynamic  analysis  in  which  the  predominant  cost  is  determined  by 
the  residual  force  calculation  and,  in  particular,  the  evaluation  of  stresses  at 
quadrature  points.  In  this  case,  storage  and  compute  cost  are  directly  pro¬ 
portional  to  the  number  of  quadrature  points.  Typical  commercial  explicit 
codes,  which  represent  the  dominant  technologies  utilized  in  crash  analysis 
and  metal  forming,  typically  employ  low-order  elements,  usually  four-node 
quadrilateral  shell  elements  with  one  through-thickness  stack  of  quadrature 
points,  and  eight-node  hexahedral  elements  in  three-dimensional  solid  anal¬ 
ysis  with  one  quadrature  point.  The  location  of  the  quadrature  points  is  at 
the  origin  of  parametric  coordinates  within  the  element,  that  is,  the  one-point 
Gauss  rule.  This  minimizes  storage  of  stresses  and  the  number  of  constitutive 
evaluations  and  results  in  an  efficient  computational  procedure  for  very  large 
industrial  problems.  However,  there  are  shortcomings  engendered  by  one-point 
quadrature  for  which  no  completely  satisfactory  solutions  have  been  found. 

The  most  prominent  shortcoming  is  that  one-point  quadrature  results  in  rank 
deficiency  of  the  discrete  system.  The  poster  child  for  this  phenomenon  is 
the  famous  “hourglass  mode” ,  although  other  modes  of  deformation  also  give 
rise  to  singularities  (see,  e.g.,  Chapters  4.6  and  5.3  in  [25]).  For  simplicity  of 
vocabulary,  henceforth  we  will  simply  refer  to  all  spurious  singular  modes  as 
hourglass  modes. 

There  has  been  a  very  large  number  of  papers  devoted  to  stabilizing  hourglass 
modes  by  artificial  viscous  and  elastic  mechanisms,  but  the  commonly  used 
techniques  are  ad  hoc  and  involve  parameters  that  require  tuning.  If  the  pa¬ 
rameters  are  too  small,  mesh  instabilities  appear,  if  they  are  too  large,  physical 
response  is  altered.  In  many  case  the  happy  medium  is  difficult  to  predict.  If 
runs  need  to  be  repeated  and  parameters  adjusted  in  trial  and  error  fashion, 
engineering  and  computer  time  are  wasted. 

Another  shortcoming  is  the  low-order  accuracy  of  standard  explicit  dynamics 
algorithms.  The  hope  of  achieving  higher-order  accuracy  in  explicit  finite  cle¬ 
ment  technology  seems  to  have  been  lost  long  ago.  Why  has  this  occurred?  In 
our  view  it  has  to  do  with  the  discrete  spectrum  of  higher-order  finite  elements. 
As  the  degree  of  polynomials  is  increased,  the  lower  modes  become  more  ac¬ 
curate,  which  is  theoretically  well  understood,  but  the  higher  modes  diverge 
with  polynomial  degree,  a  surprising  result  that  was  discovered  recently  in 
the  investigations  of  [15,28,36];  see  also  [14].  These  many  “bad”  modes  are 
not  much  of  an  issue  in  elliptic  (e.g.,  static  elastic)  or  parabolic  (e.g.,  heat 
conduction)  problems,  due  to  the  strong  stability  of  these  operators,  but  they 
are  a  significant  problem  in  hyperbolic  (e.g.,  structural  dynamic)  problems 
and,  in  particular,  ones  in  which  high  modal  behavior  is  unavoidable,  such  as 
during  impact  in  crash  dynamics.  Higher-order  finite  elements  exhibit  a  lack  of 
robustness  in  these  cases.  Low-order  finite  elements  have  the  advantage  that 


2 


their  higher  modes  are  better  behaved  than  those  of  high-order  finite  elements, 
and  this  seems  to  be  an  important  reason  for  the  preference  of  low-order  el¬ 
ements.  Poor  practical  experiences,  which  years  ago  led  to  the  elimination 
of  higher-order  elements  in  explicit  codes,  such  as  LS  DYNA,  may  also  be 
attributed  to  these  negative  results  for  higher-order  finite  elements.  However, 
based  on  recent  innovations,  we  believe  that  there  may  be  alternatives  to  what 
currently  exists. 

One  direction  that  may  be  pursued  in  an  effort  to  develop  better  minimal- 
quadrature-point  discretization  technology  is  collocation,  but  within  the  frame¬ 
work  of  Isogeometric  Analysis  (IGA),  which  has  been  the  subject  of  numerous 
recent  studies  [1-16,19,21,22,24,26,28,29,36,39,40]  ,  rather  than  within  tradi¬ 
tional  finite  element  analysis.  IGA  utilizes  smooth  basis  functions  emanating 
from  computer  aided  geometric  design,  for  example,  NURBS,  T-splines,  etc. 
The  raison  d’etre  of  IGA  is  to  simplify  the  generation  of  finite  element  mod¬ 
els  from  CAD  designs  by  utilizing  a  single  mathematical  representation  for 
both  design  and  analysis.  However,  IGA  has  also  been  shown  to  be  a  superior 
computational  mechanics  technology  to  traditional  finite  elements  in  many 
situations  [ibid.].  In  particular,  and  of  importance  in  hyperbolic  cases,  is  that, 
unlike  traditional  finite  elements,  the  higher  modes  of  IGA  basis  functions 
do  not  diverge  with  increasing  degree,  but  in  fact  achieve  almost  spectral 
accuracy  that  improves  with  degree.  It  has  been  shown  that  the  robustness 
of  higher-order  NURBS  elements  increases  with  polynomial  degree  [32],  in 
contrast  with  the  behavior  of  higher-order  finite  elements.  This  deficiency  of 
traditional  higher-order  finite  elements  is  overcome  by  IGA. 

The  smoothness  of  IGA  basis  functions  enables  use  of  the  strong  form  of  the 
partial  differential  equations,  which  provides  a  platform  for  collocation  with 
interesting  stability  properties.  Collocation  may  be  viewed  as  a  variant  of 
one-point  quadrature.  It  is  simple  to  show  that  for  quadratic  and  higher-order 
NURBS,  with  uniform  knot  vectors  and  a  suitable  choice  of  the  collocation 
points,  the  discrete  Laplace  operator  produced  by  collocation  is  rank  sufficient 
in  all  dimensions.  It  follows  that  the  elasticity  operator  is  also  rank  sufficient, 
that  is,  no  hourglass  modes.  This  is  in  contrast  with  one-point  quadrature 
on  low-order  finite  elements,  as  described  above.  The  upshot  of  this  observa¬ 
tion  is  that  IGA  collocation  methods  eliminate  the  need  for  ad  hoc  hourglass 
stabilization  techniques  and  their  tuning  parameters. 

Based  on  the  previous  discussion,  IGA  collocation  opens  the  way  to  stable,  ro¬ 
bust,  higher-order  accurate  methods  with  a  minimum  number  of  quadrature 
points.  In  order  to  take  advantage  of  these  attributes  in  explicit  dynamics, 
one  must  use  lumped  mass  matrices.  A  unique,  positive,  lumped  mass  matrix 
may  be  computed  for  IGA  by  the  row-sum  technique  (see  [25],  Chapter  7.3). 
Unfortunately,  the  lumped  mass  matrix  does  not  maintain  the  accuracy  of  the 
consistent  mass  matrix.  This  has  been  noted  in  our  previous  studies  [15,28,14], 
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In  fact,  accuracy  is  limited  to  second-order  no  matter  the  polynomial  degree, 
and  accuracy  degrades  very  significantly  in  higher  modes  (see  also  [39]).  A 
solution  to  this  dilemma  is  provided  by  the  explicit  predictor  multi-corrector 
algorithms  described  in  (see  [27],  Chapter  9.4  in  [25]  and  Chapter  6.2.3  in  [14]) 
in  which  the  lumped  mass  is  used  as  the  diagonal  “left-hand-side”  matrix  and 
the  consistent  mass  is  used  in  the  “right-hand-side”  residual  vector  calcula¬ 
tion.  This  preserves  the  usual  explicit  computational  architecture  and  through 
corrector  iterations  is  able  to  maintain  the  spatial  accuracy  of  consistent  mass. 

The  combination  of  the  technologies  described  above  provides  the  potential 
for  efficient,  stable,  robust  and  higher-order  accurate  explicit  dynamics  meth¬ 
ods.  This  paper  is  devoted  to  pursuing  the  issues  arising  in  the  development 
of  this  methodology.  It  builds  on  our  earlier  work  [3]  in  which  we  introduced 
a  collocation  scheme  for  a  single  NURBS  patch.  In  [3]  we  provided  a  com¬ 
plete  mathematical  analysis  of  the  one-dimensional  case.  The  results  of  the 
mathematical  analysis  are  not  applicable  in  multiple  dimensions  and  thus  this 
remains  an  open  problem.  However,  we  confirmed  the  theoretical  convergence 
rates  numerically  on  linear  elliptic  problems  in  one,  two  and  three  dimensions. 
We  also  numerically  investigated  the  discrete  eigenspectrum  and  the  effects  of 
alternative  locations  of  collocation  points. 

An  outline  of  this  paper  follows: 

In  Section  2  we  briefly  review  some  basic  results  on  B-splines,  NURBS  and 
the  resulting  element  structure  of  multi-patch  configurations. 

In  Section  3  we  start  with  a  standard  variational  formulation  for  the  linear 
elastostatic  problem  and,  invoking  sufficient  smoothness  of  the  function  spaces, 
integrate  by  parts  to  yield  the  strong  form  of  the  residual  within  patches  and 
traction  continuity  conditions  on  patch  interfaces  and  external  boundaries. 
We  do  not  assume  the  test  function  space  is  the  same  as  the  solution  space. 
We  then  use  standard  techniques  to  construct  C°°  test  functions  with  com¬ 
pact  support  that  comprise  Dirac  delta  sequences  about  each  of  the  desired 
collocation  points.  Taking  the  limit  of  the  sequences,  defines  the  collocation 
scheme.  Special  consideration  needs  to  be  given  to  points  on  the  patch  inter¬ 
faces  and  external  boundaries  and  we  give  precise  descriptions  of  the  treatment 
of  several  important  cases. 

In  Section  4  we  generalize  to  dynamics  and  describe  the  explicit  predictor 
multi-corrector  time  integration  algorithms.  We  argue  that  if  a  sufficient  num¬ 
ber  of  explicit  multi-corrector  iterations  is  utilized,  the  higher-order  spatial  ac¬ 
curacy  of  the  corresponding  implicit  algorithm  with  consistent  mass  is  achieved. 

In  Section  5  we  present  several  static  numerical  examples.  We  test  the  method 
on  Dirichlet  and  Neumann  boundary  conditions,  mixed  boundary  conditions, 
and  on  single  and  multi-patch  configurations  with  different  material  properties 
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in  the  patches.  Satisfactory  results  are  obtained  in  all  cases. 

In  Section  6  we  present  dynamic  cases.  We  confirm  the  higher-order  con¬ 
vergence  rates  of  the  explicit  multi-corrector  method  on  a  one-dimensional 
example  and  a  two  dimensional  plane  strain  annular  configuration. 

We  draw  conclusions  in  Section  7. 


2  NURBS-based  isogeometric  analysis 


Non-Uniform  Rational  B-splincs  (NURBS)  are  a  standard  tool  for  describing 
and  modeling  curves  and  surfaces  in  computer  aided  design  and  computer 
graphics  (see  Piegl  and  Tiller  [34]  and  Rogers  [37]  for  an  extensive  description 
of  these  functions  and  their  properties).  In  this  work,  we  use  NURBS  as  an 
analysis  tool,  as  proposed  by  Hughes  et  al.  [26].  The  aim  of  this  section  is 
to  present  a  short  description  of  B-splines  and  NURBS,  followed  by  a  simple 
discussion  on  the  basics  of  isogeometric  analysis  and  by  an  introduction  to  the 
proposed  collocation  method. 


2.1  B-splines  and  NURBS 


B-splines  in  the  plane  are  piecewise  polynomial  curves  composed  of  linear 
combinations  of  B-splinc  basis  functions.  The  coefficients  (B<:)  are  points  in 
the  plane,  referred  to  as  control  points. 

A  knot  vector  is  a  set  of  non- decreasing  real  numbers  representing  coordinates 
in  the  parametric  space  of  the  curve 

{6  =  0,...,^+p+i  =  l},  (1) 

where  p  is  the  order  of  the  B-spline  and  n  is  the  number  of  basis  functions 
(and  control  points)  necessary  to  describe  it.  The  interval  [£i,  £n+p+i]  is  called 
a  patch.  A  knot  vector  is  said  to  be  uniform,  if  its  knots  are  uniformly-spaced 
and  non-uniform  otherwise;  it  is  said  to  be  open  if  its  first  and  last  knots 
have  multiplicity  p+  1.  In  what  follows,  we  always  employ  open  knot  vectors. 
Basis  functions  formed  from  open  knot  vectors  are  interpolatory  at  the  ends  of 
the  parametric  interval  [0, 1]  but  are  not,  in  general,  interpolatory  at  interior 
knots. 

Given  a  knot  vector,  univariate  B-spline  basis  functions  are  defined  recursively 
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starting  with  p  =  0  (piecewise  constants) 


i  if  6  <  £  <  6+i 

0  otherwise 


(2) 


For  p  >  1  : 


Ni,P{0  =  {  6  •  P  -  6 

0 


^  +  /i+P+1  /  Ni+1,p_ x(0 


6+p+i  6+i 


where,  in  (3),  we  adopt  the  convention  0/0  =  0. 


if  6  <  £  <  6+p+i 
otherwise, 

(3) 


In  Figure  1  we  present  an  example  consisting  of  n  =  9  cubic  basis  functions 
generated  from  the  open  knot  vector  {0,  0, 0,  0, 1/6, 1/3, 1/2,  2/3,  5/6, 1, 1, 1, 1}. 


Fig.  1.  Cubic  basis  functions  formed  from  the  open  knot  vector 
{0, 0, 0, 0, 1/6, 1/3, 1/2, 2/3, 5/6, 1, 1, 1, 1}. 


If  internal  knots  are  not  repeated,  B-spline  basis  functions  are  Cp~l -continuous. 
If  a  knot  has  multiplicity  k,  the  basis  is  Cp_fc-continuous  at  that  knot.  In  par¬ 
ticular,  when  a  knot  has  multiplicity  p,  the  basis  is  C°  and  interpolates  the 
control  point  at  that  location.  We  define 

S  =  span{W,P(06  =  1,  •  •  ■ ,  n}  (4) 


By  means  of  tensor  products,  a  multi-dimensional  B-spline  region  can  be  con¬ 
structed.  We  discuss  here  the  case  of  a  two-dimensional  region,  the  higher¬ 
dimensional  case  being  analogous.  Consider  the  knot  vectors  {£i  =  0, ...,  6i+P+i 
1}  and  {?7i  =  0,  ...,rjm+q+i  =  1},  and  an  n  x  rn  net  of  control  points  B,j.  One 
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dimensional  basis  functions  Ni)P  and  Mhq  (with  i  =  1, n  and  j  =  1, m)  of 
order  p  and  q,  respectively,  are  defined  from  the  knot  vectors,  and  the  B-spline 
region  is  the  image  of  the  map  S  :  [0, 1]  x  [0, 1]  — >•  Q  given  by 

n  m 

s«,  l)  =  EE  WmH  By.  (5) 

i= 1  3= 1 

The  two-dimensional  parametric  space  is  the  domain  [0,1]  x  [0,1].  Observe 
that  the  two  knot  vectors  {£i  =  0, ...,  £n+p+i  =  1}  and  {rji  =  0,  ...,r)m+q+i  =  1} 
generate  a  mesh  of  rectangular  elements  in  the  parametric  space  in  a  natural 
way.  Analogous  to  (4),  we  define 

S  =  span {Nij,(£)Mjiq(r}),i  =  1, . . . ,  n,  j  =  1, . . . ,  m}  (6) 


In  general,  a  rational  B-spline  in  is  the  projection  onto  d-dimensional  phys¬ 
ical  space  of  a  polynomial  B-spline  defined  in  (d+l)-dimensional  homogeneous 
coordinate  space.  For  a  complete  discussion  see  the  book  by  Farin  [23]  and 
references  therein.  In  this  way,  a  great  variety  of  geometrical  entities  can  be 
constructed  and,  in  particular,  all  conic  sections  in  physical  space  can  be 
obtained  exactly.  The  projective  transformation  of  a  B-spline  curve  yields  a 
rational  polynomial  curve.  Note  that  when  we  refer  to  the  “degree”  or  “or¬ 
der”  of  a  NURBS  curve,  we  mean  the  degree  or  order,  respectively,  of  the 
polynomial  curve  from  which  the  rational  curve  was  generated. 

To  obtain  a  NURBS  curve  in  M2,  we  start  from  a  set  e  M3  (i  =  1,  ...,n) 
of  control  points  (“projective  points”)  for  a  B-splinc  curve  in  M3  with  knot 
vector  S.  Then  the  control  points  for  the  NURBS  curve  are 


Bj 


1,2 


(7) 


where  [B,]*,  is  the  kth  component  of  the  vector  B;  and  u>i  =  [B^  is  referred 
to  as  the  ith  weight.  The  NURBS  basis  functions  of  order  p  are  then  defined 
as 


m) 


(8) 


The  NURBS  curve  is  defined  by 


CK)=E«Bi. 

i=l 


(9) 


Analogously  to  B-splines,  NURBS  basis  functions  on  the  two-dimensional 
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parametric  space  fl  =  [0, 1]  x  [0, 1]  are  defined  as 


Ni,p{0MiMui,j 


*,3 


E?,i  IX, 


(10) 


where  ujh3  =  (B£  -  )3.  Observe  that  the  continuity  and  support  of  NURBS  basis 
functions  are  the  same  as  for  B-splincs.  NURBS  spaces  are  the  span  of  the 
basis  functions  (10). 


NURBS  regions,  similarly  to  B-spline  regions,  are  defined  in  terms  of  the 
basis  functions  (10).  In  particular  a  single-patch  domain  is  a  NURBS  region 
associated  with  the  n  x  m  net  of  control  points  B4iJ,  and  we  introduce  the 
geometrical  map  F  :  D  — >  D  given  by 

n  m 

(11) 

i=l  j= 1 


In  order  to  deal  with  multi-patches  geometries,  consider  N  spline  and  NURBS 
spaces,  defined  on  U,  possibly  from  different  knot  vectors  {£1^  =  0, ...,  (nk+pk+\,k 
1}  and  {rjhk  =  0, ...,  r)mk+qk+1,k  =  1}  and  weights.  Denote  by  Efffjf  (£,  rj) 
the  basis  functions  on  the  k-th  space  and  consider  then  the  NURBS  maps 
F  ^  — y  klk 


Tlk  '171  k 

F*~E£  T((,v)V,d,t,  k  =  l,2,...,N, 

i= 1 j= 1 

where  the  Qk  C  M2  represent  the  patches  which  compose  the  connected 
(closed)  physical  domain 

U  =  (J  Qk. 

k=l,2,...,N 

By  D  we  indicate  the  interior  part  of  D.  We  assume  that  the  maps  F*,  e  C^fl) 
and  their  inverses  F^1  G  C,1(h2fc).  Moreover  we  assume  that  the  patches  Qk 
form  a  geometrically  conforming  multi-patch  structure;  see  for  instance  figure 
3. 

Assumption  2.1  (Geometrical  conformity)  The  common  boundary  of  any 
two  patches  is  either  void,  or  the  image  of  corners,  or  the  image  of  edges  of 
12.  More  precisely,  we  have 

dVLk  n  dVLk,  =  F k(E)  =  F kfE')  1  <  k,k'  <  N, 

with  E,  E 7  either  two  (possibly  coincident)  corners  or  two  (possibly  coincident) 
edges  of  fl  . 


Fig.  2.  Example  of  a  domain  described  by  three  patches. 


Remark  2.1  The  extension  to  three  dimensions  is  analogous  and  not  dis¬ 
cussed. 

2.2  NURBS  on  the  physical  domain 

Following  the  isoparametric  approach,  the  space  of  NURBS  vector  fields  on 
each  patch  iR  is  defined,  component  by  component  as  the  span  of  the  push- 
forward  of  the  basis  functions  (10) 


We  assume  that  the  spaces  are  conforming  at  the  patch  boundaries. 

Assumption  2.2  (Space  conformity)  With  the  notation  of  Assumption  2.1, 
the  spaces  V h,k\Fk(E )  and  Nh,k'\'Fk,(E')  coincide. 

Then,  by  continuous  gluing  of  (12),  we  define 


(13) 


The  image  of  the  elements  in  the  parametric  space  are  elements  in  the  physical 
space.  The  elements  of  the  physical  mesh  in  are  therefore 


%  =  {Ffc((&,£i+i)  x  with  i  =  l,...,nk+pt,  j  =  1,.. .  ,mk  +  • 


(14) 


We  denote  by  7  =  14=1,2 and  by  h  the  mesh-size,  that  is,  the  maximum 
diameter  of  the  elements  of  7. 


9 


The  interested  reader  may  find  more  details  on  isogeometric  analysis  as  well  as 
many  interesting  applications  in  a  number  of  recently  published  papers  [6,15,16,14,26]. 

Remark  2.2  Let  E  and  E'  be  two  edges  as  in  Assumption  2.1.  Then,  as  a 
consequence  of  Assumption  2.2,  the  elements  generated  on  the  common  bound¬ 
ary  d Llk  H  dLlk1  by  the  mappings  F k  and  F k>  coincide.  The  same  holds  for  the 
two  associated  knot  vectors. 


3  Isogeometric  collocation  method  for  elastostatics  in  a  variational 
context 


In  this  section  we  follow  a  variational  construction  in  order  to  derive  the  col¬ 
location  equations.  In  order  to  fix  ideas,  we  will  concentrate  our  presentation 
on  the  linear  elasticity  problem. 


3.1  Problem  description  and  variational  formulation 


Let  C  represent  an  elastic  body  that  is  subjected  to  prescribed  displace¬ 
ment  on  part  of  its  boundary  T n,  and  to  (possibly  null)  traction  h  G  [L°°(r 7v)]d 
on  the  remaining  part  of  the  boundary  T n.  We  assume  that  T n  and  T n  are 
made  of  a  finite  union  of  connected  and  regular  components.  We  assume  that 
the  traction  h  is  piecewise  continuous,  i.e. ,  is  continuous  on  the  image  of  each 
face  (or  edge  if  d  —  2)  of  O.  Regarding  the  boundary  conditions  on  T^,  we 
make  use  of  a  function  g  G  [C°(rjD)]d  as  shown  below.  The  body  is  also 
subjected  to  a  volume  loading  f  G  [L°°(f2)]d  such  that  f|nfc  G  [C'°(f2/c)]d  for 
k  =  1,2,  ...,1V.  We  indicate  with  C  the  standard  fourth-order  elasticity  ten¬ 
sor.  We  assume  that  all  components  of  C*,  :=  C|nfc  are  in  C'1(f4),  while  we 
allow  C  to  have  jumps  from  one  patch  to  another.  The  above  regularity 
requirements  on  the  data  can  be  made  weaker;  see  Remark  3.1  below. 

Given  the  function  spaces 

Vg  =  {„  6  [tf ‘(Si)]*  :  t>|r„  =  g}, 

V„  =  [v  6  =  0}, 

the  elasticity  problem  in  variational  form  (based  on  the  principle  of  virtual 
work)  reads 

{Find  u  G  14  such  that 

f  CVs(u)  :  V5(v)  =  f  f  •  v  +  /  h  •  v  Vv  G  Vq. 

J  J  £1  ifiu 
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3.2  A  family  of  discrete  problems 


In  order  to  simplify  the  exposition,  we  extend  the  function  h  by  zero  on  the 
whole  Ufc=i ,2,...,iv<9^fc-  Thus  h  represents  a  vector  function  living  on  Ufc=i 
whose  support  is  contained  in  T jy.  We  now  assume  that  the  solution  u  is  in 
C2(flfc)  for  all  k  =  1,2,  Then,  integrating  by  parts  on  each  patch  Qk, 

k  =  1,2,  ..,N,  and  rearranging  terms,  equation  (15)  gives 

I  (div  CfcV^u)  +  f)  •  v  +  ^  [  (cfeV5'(u)  ■  nfc  -  hV  v  =  0  ,  (16) 
Jnl  J  Jaul  J 

for  all  v  G  Vo,  where  n^.  represents  the  outward  unit  normal  to  the  domain 

klk- 

We  will  now  consider  for  simplicity  the  two-dimensional  case  d  —  2,  while  the 
extension  to  three  dimensional  problems  will  be  discussed  afterwords. 

We  consider  the  case  where  the  collocation  points  are  chosen  using  the  (tensor- 
product)  Greville  abscissae  [18,20]  of  the  knot  vectors.  What  follows  applies 
also  to  other  sets  of  collocation  points  such  as,  for  instance,  the  Demko  ab¬ 
scissae,  see  [20]. 

Let  k  G  {1,  2,  Let  k,  i  —  1, ...,  rik,  be  the  Greville  abscissae  related  to 

the  knot  vector  ...,£nk+Pk+1  ifc}: 

£i,fc  =  (£i+l,fc  +  £i+2,k  +  •••  +  fi+pk,k)/Pk  ■  (17) 

Analogously,  let  j  =  1  be  the  Greville  abscissae  related  to  the 

knot  vector  {qlik,  ...,pmk+qk+ltk}.  It  is  easy  to  see  that  f  l  k  =  rjl  k  =  0,  lnk  k  = 
Vmk  k  =  I;  while  all  the  remaining  points  are  in  (0, 1).  We  define  the  collocation 
points  %,j,k  G  O  by  the  tensor  product  structure 

A.j./c  \^i,ki  Vj,k)  G  kl  , 

for  i  =  1 ,...,nk,  j  =  1 , . . . ,  mk ■  The  collocation  points  on  12  are  then  defined 
accordingly  for  all  k  =  1, 2, ...,  N  and  all  i ,  j  as  above 

F k(d~i,j,k)  G  12fc  C  12  . 

Note  that  collocation  points  on  <9f 4  fl  dQj..  ^  0  coincide  because  of  our  con¬ 
formity  assumptions,  see  Remark  2.2. 

Let  functions  5f^k  G  G°°(r2),  for  e  G  M+  and  for  the  same  sets  of  indices  i,j,  k, 
as  given  above,  be  defined  as  follows.  Let  p  :  M+  — y  M  be  the  C°°  positive 
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function 


ip(x)  = 


e±/{x2-i)  x  e  [0)  i) 

x  o  x  g  [l,  +oo)  . 

We  then  define  the  C°°  radial  functions  :  R2  — y  M+,  depending  on  the  real 

parameter  0  <  £  <  1,  as 


^£(x)  =  p(||x||/e)  VxgM2, 

where  ||  ■  ||  indicates  the  euclidean  norm.  Finally,  given  any  collocation  point 
Tijtk  G  D ,  for  k  =  1,2,...,  At  and  i  =  1  ,...,nk,  j  =  1  ,...,rrik,  we  dehne  the 
functions  8£l  ]  k  as 


^j,fe(x)  =  ^£(x  “  T<J,fc)  (  Jn  VE  -  TiJ,k) 


-i 


VxeO. 


(18) 


As  is  evident,  d|j  fc  is  a  smoothed  Dirac  delta  function  converging  to  the  Dirac 
delta  distribution  located  at  as  £  — *  0. 

From  Lemma  A.l  it  follows,  for  all  s  G  [L°°(D)]2  that  are  continuous  about 
each 

s  5tj,k  =  s(Ti,j,k ),  (19) 

where  the  integral  above  is  to  be  calculated  component  by  component. 


In  order  to  derive  a  collocation  method  from  the  variational  equation  (15),  we 
search  for  a  bounded  function 


U/i  G  g  {v/j  G  V/j  .  V hij~i,j,k)  ^^~i,j,k  G 

that  satisfies  (16)  for  all  vector  test  functions  [0,  <5|Gfc]  and  [5|j  fc,0],  with  k  = 
1,2,  ...,1V  and  i  =  l,...,n*,,  j  =  1,  ...,mfc.  We  assume  that  Vh  C  C°(fl)  D 
C2(Ufcf2fc),  where  by  C2(Uk^ik)  we  mean  that  u^|nfc  is  of  class  C2,  for  each 
k  =  1,2,...,  N. 

Note  that,  due  to  the  definition  of  the  discrete  space  Vh,  we  are  strongly 
enforcing  the  global  C°  continuity  of  the  solution  in  fl. 


Testing  on  the  above  functions  gives  the  following  set  of  equations  depending 
on  the  parameter  £ 


N  r  ( 

-E  (dwC-kWs(nh) 

k=iJsi* 


f )  + 


N 

E 

k=  1 ' 


'on 


aV5(u h)  ■  -  h 


=  o 


(20) 

for  k  —  1,2,  i  =  1, ... ,nk ,  j  =  1,  and  where  as  usual  an  integral 

of  vector  functions  is  calculated  component  by  component. 
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3.3  Collocation  equations 


The  collocation  method  is  obtained  in  this  section  by  taking  the  limit  £  — >■  0 
of  (20).  We  first  consider  the  simpler  case  of  collocation  points  that  are  in  the 
interior  of  a  patch,  that  is  Tjj y  G  and  then  address  the  more  subtle  case 
of  collocation  points  that  belong  to  the  patch  boundary,  that  is  r.hJ^  G  d f^. 


3.3.1  Collocation  equations  at  the  patch  interior 

Let  Tij'k  G  As  noted  above,  such  condition  is  met  if  and  only  if  i  ^  0,  j  ^  0 
and  i  ^  nk,j  ^  .  In  such  case,  since  the  support  of  Sfj  k  is  contained  in 

the  ball  of  radius  £  centered  at  t*^,  the  second  term  in  (20)  vanishes  for 
sufficiently  small  e.  Moreover,  for  the  same  range  of  £,  all  the  terms  but  one 
in  the  first  sum  of  (20)  vanish,  thus  giving 

/n(divC»vs(uO  +  f).yM  =  o. 

Taking  the  limit  for  £  — >■  0  in  the  above  equation  and  recalling  (19),  we  get 

^div  CfcVs(u/l)  +  f^)  (Tjj, k)  =  0,  (21) 

that  is,  the  collocation  of  the  strong  form  of  the  equations  at  Tjjy. 


3.3.2  Collocation  equations  at  the  patch  boundary 

Assume  now  rhhk  G  d Vtk-  We  follow  the  rules: 

•  if  Ti,j,k  ^  T £)  we  do  not  set  a  collocation  equation  (however,  we  have 
u h{ji,j,k)  =  g (Ti,j,k),  since  the  Dirichlet  boundary  condition  is  enforced  a 
priori  in  W)g); 

•  each  Tijtk  G  Ttv  is  associated  with  a  collocation  equation,  that  sets  the 
traction  on  the  boundary  T N ; 

•  the  remaining  Tjjy  G  <9f4  belong  to  the  inter-patch  boundary,  and,  as 
already  observed;  these  coincide  with  other  boundary  points  from  other 
patches;  introducing 

X  =  {kN  :  l  <k  <  N,  Tij'k  G  dQj.}  .  (22) 

Then  #%,  the  cardinality  of  X,  is  the  number  of  patches  that  share  | 
in  this  situation  we  associate  to  Tjjy  a  single  collocation  equation  as  shown 
below.  Note  moreover  that,  since  the  space  g  is  continuous  a  priori,  there 
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holds  implicitly  —  1  gluing  conditions,  namely 


ufc|nfel  (Ti,j,k)  —  u/i|ofc2  iTi,j,k )  —  •  •  •  —  (Tj^fc),  for  ki,  k2,  ■  ■  ■ ,  G  X 


Considering  Section  3.3.1  and  the  rules  above,  it  is  easy  to  check  that  the  final 
number  of  collocation  equations  corresponds  to  the  size  of  the  discrete  space 


Considering  then  Ttj,k  G  dflk\^D,  we  now  derive  the  collocation  equation, 
passing  to  the  limit  in  (20).  Many  possibilities  arise,  some  of  them  are 

1)  point  Tijtk  is  on  the  boundary  of  a  single  patch,  G  dQ  and  it  is  not 
a  corner; 

11)  point  Tijtk  is  on  the  boundary  of  a  single  patch  and  is  a  corner  of  fl; 

III)  point  Tij}k  belongs  to  the  common  boundary  of  two  patches  Qk  and  is  not 
a  corner  of  both  of  them  (and  therefore  is  inside  Q); 

IV)  point  Tjjy  is  a  common  corner  of  two  or  more  patches  f 1*,; 

and  are  represented  in  Figure  3. 


O  Dirichlet  81  Type  III 

•  Type  I  □  Type  IV  (boundary) 

I°1  Type  II  ■  Type  IV  (internal) 


Fig.  3.  Example  of  collocation  point  types  according  to  Table  1. 


The  collocation  equations  in  these  cases  are  reported  in  Table  1,  and  are 
derived  by  the  following  general  approach.  Given  Tjjy,  %  as  defined  in  (22) 
and  any  k  G  X,  for  £  sufficiently  small  the  curve  dflk  FI  supp(<5|Ffc)  can  be 
divided  into  two  C^-connected  components  E 'k  and  E”,  such  that 


dn~k  n  supp  (Sljk)  =  E'kUE] 


i// 

k  > 


(23) 
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Case 

Equation 

I)  Ti,j,k  G  not  a  corner  point 

(CfcVs(u/,)  ■  nk  -  h)  =  0 

II)  Tijfi  G  a  corner  point,  %  =  {k} 

(CfcV^fu/,)  •  n'k  -  h')  (u j,fc)  +  (CfcV^fu/,)  ■  n'l  -  h")  =  0 

III)  Tijtk  e  dflk  n  dQki  Tj  £  dfl,  X  =  {k,  k} 

(CfcVRuhlnJ  ■  n k)(Ti,j,k)  +  (C^.Vs(uh|nE)  •  n k)(n,j,k)  a  [cVs(uft)  ■  n|  (Tij,fc)  =  0 

IV)  on  the  boundary: 

Ti,j,k  G  dQk  n  <9%,  Tijtk  €  <9ft,  X  =  {k,k} 

(CfcVs(uh)  ■  n'j,  -  h)  (7ij,fc)  +  (ck Vs (uh)  ■  nC  -  h)  (n^k) 

+  |cVs(uh)  ■  nj  (Tijtk)  =  0 

IV)  internal  to  the  domain: 

ri,j,k  £  p|  a  corner  point, 

eex 

Ti,j,k  i  dfi,  X  =  {k,  k,  k} 

Z^e3c(cWS(uh|nd  •  rie)(Tijjtk)  +  (QVs(uh|nd  ■  n ")(nj,fc)  =  0 

Table  1 


Some  examples  of  the  boundary  equations,  i.e.,  assuming  ryyfc  £  dQk  with  $ 
To-  We  refer  to  Figures  3  and  4  for  a  graphical  representation  of  the  tabulated 
examples.  The  symbol  [  ]  stands  for  the  standard  jump  operator  with  n  indicating 
a  unit  normal  to  dflk  n  The  first  and  second  cases  are  examples  of  Neumann 
boundary  conditions.  The  third  case  represents  a  gluing  condition  between  two 
patches  (here  Tijk  is  not  a  corner),  while  the  fourth  and  fifth  are  a  mix  of  the 
above.  A  description  of  the  different  cases  is  presented  in  the  text. 

Note  that  for  simplicity  of  exposition  we  have  left  implicit  the  dependence  of 
Et,  E'i  on  e.  Multiplying  (20)  by  the  positive  parameter  e  yields,  for  sufficiently 
small  e, 

-£  £  /  (div  CiVs(uh)  +  f )  +  £  E  /  (CiVS(<M)  •  ns  -  h)  6!j_t  =  0  . 

Ser™^  '  SeJC  Jaa*  ' 

(24) 

Due  to  Lemma  A.l  m  the  appendix,  for  all  k  e  %  the  integral  /„l(divClV>J1)+ 

f  j  Sfj  k  is  bounded  independently  of  e.  Therefore  the  first  term  in  (24)  vanishes 

in  the  limit  e  — >  0.  Regarding  the  second  term  in  (24),  recalling  (23)  for  all 
A;  £  IK  we  can  write 

^  /  (CjfcV5(uft)  ■  n£  -  h)  S;jk 

= £  X(CsvS(ufe)  ■  n~k  - h)  s^k + £X(c^v5(u/i)  ■ ns  _  h)  ^  ■ 

k  k 

(25) 

Note  that,  if  Tlj  k  is  a  corner  point,  the  two  (outward)  normal  vectors  at  Tl)])k 
are  different  for  E'r  and  E'l.  We  therefore  indicate  with  nr  the  unit  outward 
(with  respect  to  ttk)  normal  to  E'k  calculated  in  Tk]^-  Analogously,  we  indicate 
with  n g  the  unit  outward  normal  to  Ej'  calculated  in  (see  Figure  4). 
Similarly,  since  the  datum  h  is  allowed  to  jump  at  we  define  h'  and  h" 
as  the  limit  values  of  h  at  obtained  respectively  for  E'r  and  E'-l. 
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Taking  the  limit  and  using  Lemma  A. 2  thus  yields 


tone  JJc-kX7s(uh)  ■  n'k  -  h')  Sfdtk  =  Cs  (ckVs(uh)  •  n'k  -  h')  (riJik)  ,  (26) 

k 

and  the  analogous  result  holds  for  E'l.  Combining  all  the  previous  arguments, 
taking  the  limit  in  (24),  using  (25),  (26)  and  dividing  the  result  by  Cs,  finally 
gives 


E 


n'fc  -  h/)(Ay,fc)  +  ( CkVs(uh ) 


0  .  (27) 


Equation  (27)  enforces  the  Neumann  boundary  condition  on  T jy  and  the  nor¬ 
mal  stress  continuity  at  the  inter-patches  boundaries. 


3.3.3  Final  system  of  equations  at  the  collocation  points 

The  summary  of  the  equations  at  the  collocation  points  Tyy  is  reported  in 
Table  2. 

The  final  linear  system  reads  as  usual 

Ku  =  F,  (28) 

where  K  denotes  the  global  stiffness  matrix,  F  incorporates  the  body  force 
(f)  and  inhomogeneous  Diri chief  (g)  and  Neumann  (h)  boundary  conditions 
(see  [25]),  and  u  is  the  unknown  vector  of  control  variables  associated  to  the 
discrete  displacement  (u/J.  It  is  convenient  to  compute  (28)  by  assembling 
the  contributions  from  the  various  patches.  Indeed,  for  all  Tjjy  G  dQkl  k  = 
1,  2,  ...,N,  equation  (28)  is  nothing  but  the  sum  of  the  contributions  from  each 
different  patch.  This  implementation  aspect  is  also  discussed  in  Section  5.3. 


Position 

Equation 

^  ^ k 

(div  CfcVAu/,)  +  f)  (Tjjy)  =  0 

Efcenc  (Cs Vs(uft)  ■  nt  -  h')  (r^y)  +  (ckX7s(uh)  •  -  h")  (njy)  =  0 

Ti,j,k  ^  r D 

U h{n,j,k)  =  g (n,j,k),  included  in  the  def.  of  Vyg 

Table  2 


Equations  at  the  collocation  points  r,jy,  for  k  =  1,2,...,  A  and  i  =  1  ,...,nk,  j  = 
1,  ...,mk.  The  set  %  is  defined  in  (22).  We  recall  that  h  =  0  in  f2\rjv,  and  that  for 
coinciding  points  (ryy  =  Tj jy,  k  /  k)  the  collocation  equation  is  included  in  the 
global  system  only  once. 

Remark  3.1  In  order  to  apply  the  above  arguments,  we  assumed  that  the 
data  f ,  h  are  piecewise  C°  and  that  C  is  piecewise  C1  and  globally  C°  on  Q. 
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121  Type  II  □  Type  IV  (internal) 

Fig.  4.  Example  of  some  collocation  points  with  associated  normals. 


This  continuity  conditions  on  the  data  h,  f  and  on  C  are  introduced  only  for 
simplicity  of  exposition  and  can  be  made  weaker.  In  more  general  cases,  an 
averaged  value  of  the  problem  equation  at  the  collocation  point  needs  to  be 
considered. 


Remark  3.2  One  can  use  the  same  equations  listed  in  Table  2  in  the  three 
dimensional  case.  A  variational  justification  similar  to  the  one  we  have  pre¬ 
sented  in  2D  can  be  developed  also  for  this  case. 


4  Explicit  time  integration  for  elastodynamics 


In  this  Section  we  consider  and  discuss  the  extension  of  our  collocation  for¬ 
mulation  to  the  elastodynamic  case,  as  well  as  the  explicit  predictor  multi¬ 
corrector  algorithm  that  we  employ  for  time  integration. 
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4-1  Problem  description  and  its  algebraic  formulation 


The  elastodynamic  problem  reads,  in  variational  form, 


|  Find  u (t)  G  Vs  such  that  \/t  G  [0,T] 

I  L  CVS(UW)  :  v-v> +  L  '  (w^m)  • v = X f  ■ v + 1 h  ■ v  W  6  ^ 

(29) 

where  p  denotes  the  material  mass  density.  The  equation  is  the  analogue  of 
(15).  Discretizing  by  collocation  in  space,  as  detailed  in  Section  3  for  the  static 
problem,  we  are  led  to  the  semidiscrete  system  of  equations 


Ku  +  M 


d2  _ 
dt2U~ 


F 


where  u  =  u(f)  is  the  vector  of  displacement  control  variables  of  the  discrete 
(in  space)  displacement  uh(t,  •),  M  denotes  the  mass  or  collocation  matrix, 
K  the  stiffness  and  F  the  forcing  vector  (load  plus  boundary  conditions  of 
Dirichlet  and  Neumann  type)  in  compact  matrix  form.  Viscous  effects  may  be 
included  as  well. 


Initial  conditions 


u|t=o  —  Uq, 


d 

diu't=0 


Vo 


are  assigned.  See  [25]  for  further  details. 


4-2  Time  integration  algorithm 


We  consider  the  following  explicit  predictor  multi-corrector  time  integration 
algorithm  (see  [25]):  given  the  acceleration  aff  'pl  =  an,  velocity  v^xp*  =  vn  and 
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displacement  u^pl  =  un  at  time  t  =  nAt,  we  compute  for  r  >  1 


v„+i  =  vn  +  At(l  -  7)an, 

=  un  +  A tvn  +  (1  -  2/3)—  an, 
for  i  —  0, . , . ,  r  —  1 

MexplAa{l)  =  F  -  Ma^  -  Ku'^,  (30) 

-(*+1)  _  -(*)  ,  A  J) 

“■n+ 1  dn+ 1  W  j 

vjj i  =  vi°]i  +  yAta^J?, 

Un+i’  =  uSr  +  /5Af2a,()t+11), 

end 


where  M f,'p;  denotes  the  lumped  mass  matrix.  Assuming,  for  the  sake  of  sim¬ 
plicity,  a  constant  density  p,  due  to  the  partition  of  unity  property  of  NURBS 
the  lumped  mass  matrix  becomes  Mexpl  —  p  1,1  being  the  identity  matrix.  The 
load  vector  F  is  evaluated  at  time  t  —  (n  +  l)At.  We  refer  to  r  as  the  num¬ 
ber  of  corrector  passes.  The  algorithm  produces,  after  r  passes,  the  following 
approximations  at  time  t  —  (n  +  1)  At 


-expl  _  _ (r) 


ln+ 1 


=  a 


n+ 1) 


—expl  _  — (r) 


'  n+1 


=  V 


n+1 ) 


—expZ  _  _(r) 


U 


n+1 


=  U 


n+1- 


(31) 


From  (30)  we  derive 

aS?  =  agi  +  As«  =  (I  -  +  (M“»')_1(F  -  Kug,) 

=  (I  -  (M“>V(M  +  /SAi+JJag,  +  (M*a^*)-1(F  -  Kug,). 

Iterating,  and  using  (I  — A)_1v  =  X)[=o  A*v  +  0(||Arv||0O)  (where  ||  ■  denotes 
the  vector  maximum  norm),  and  recalling  Mexpl  —  p  I, 

=  ai+1  =  E(l  -  (M"7_1(M  +  -  Kug.) 

i= 0 

=  (M  +  /3Ai2K)-1M“p'(M“',l)-1(F  -  Kug,) 

+  0(||(l  -  +  /3A«2K))r(M“>’i)-1(F  -  Kng^Hoo) 

=  (M  +  /3A12K))-1(F-Ku«°>1) 

+  0(||(l  - />-1(M  +  /9A«2K))V‘(F  -  KugOlU). 

(33) 
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The  Newmark  algorithm  can  be  written  in  a  similar  way:  setting  a^mk  =  an, 
v%mk  =  vn  and  u^mk  =  ura,  the  time  step  is 


=  vn  +  Af(l  -  7)an 
uJJJi  =  u„  +  A tvn  +  (1  -  2f3)—an 
(M+/3At2K)Aa(0)  =  F  -  Ku^ 


Observe  that  from 


—Nmk 
An+ 1 

=  a 

—Nmk 

Vn+1 

_  ^ 

—Nmk 

— 

.  U77+l 

=  u 

(34)  we  get 

n+ 1 
(0) 


-Nmk 


n+i  +  7  At  Aan+1 


/2  a  -prNmk 


71+ 1 


+  f3At  Aa„+1 


=  (M  +  /JA«2K))-'(F  -  KuSi). 


(34) 


(35) 


In  what  follows  we  present  a  study  whose  aim  is  to  obtain  some  indication  of 
the  behavior  of  the  explicit  scheme  (30). 


It  is  well  known  that  the  Newmark  algorithm  is  second-order  accurate  in  time 
when  7  =  1/2  (see  [25])  and  as  accurate  with  respect  to  h  as  the  discretization 
scheme  in  space  allows.  Therefore,  for  the  purpose  of  investigating  the  order 
of  the  algorithm  (30),  we  compare  the  solutions  of  (30)-(31)  and  (34)-(35) 
assuming  that  u*mk  =  uenxpl  =  un,  v%mk  =  wexpl  =  vn,  a%mk  =  a(+  =  an. 
From  (33)  and  (35)  we  evaluate  a  sort  of  truncation  error 

70(||(l-p-1(M  +  /3Ai2K))’>-1(F-Ku71)||oo) 

f  (36) 

^OdlBV'tF-  KuS.JIU), 

where  the  matrix 

B  =  (I  -  p-xM)  -  p~1/3At2K  =  p-1  (  (Mexp/  -  M)  -  /3At2K^j .  (37) 

Recalling  definitions  (30),  (31)  and  (34),  we  also  get 


1 1  T^Nmk  _  ^ expl  n 
lldn+l  dn+ll|oo 

At 


YrNmk 


u; 


71+1 


— expl  1 1 
Un+1 1  loo 


At 


/3At\\a 


Nmk 

71+1 


—expl  1 1 
an+ 1 1  loo 


AtOdlBV'fF-KuSjlU), 

(38) 


In  principle,  due  to  the  structure  of  B  shown  in  (37),  we  expect  the  vector  Brw 
to  converge  to  zero  for  (h,  At)  —t  (0, 0),  at  least  when  the  vector  w  represents 
a  regular  vector  field.  Moreover,  the  larger  is  the  exponent  r,  the  higher  the 
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convergence  rate  expected.  Therefore,  for  sufficiently  large  values  of  r,  the 
explicit  method  is  expected  to  behave  like  the  Newmark  method.  Numerical 
tests  in  the  next  section  will  show  that  two  passes ,  i.e.,  r  =  2,  are  in  general 
sufficient  to  obtain  an  optimal  behavior  of  the  error  in  h,  at  least  for  the  cases 
when  p  <  5. 


5  Numerical  tests:  statics 


In  this  section  we  present  the  numerical  results  relative  to  several  2D  plane- 
strain  elasto-static  problems,  all  exploiting  the  proposed  collocation  approach, 
addressing  different  aspects  such  as  Dirichlet  versus  Neumann  boundary  con¬ 
ditions,  single-patch  versus  multi-patch  problems,  etc.  The  investigated  prob¬ 
lems  are  organized  as  follows: 

•  clamped  quarter  of  an  annulus,  testing  a  single-patch  solution  and  Dirichlet 
boundary  conditions; 

•  traction  test  for  a  single  material,  testing  a  single-patch  solution  and  mixed 
boundary  conditions; 

•  traction  test  for  a  composite  material,  testing  a  multi-patch  solution  and 
mixed  boundary  conditions. 

In  the  following  we  discuss  in  details  each  problem. 


5. 1  Clamped  quarter  of  an  annulus 


We  consider  a  quarter  of  an  annulus,  as  sketched  in  Figure  5,  with  internal 
and  external  radii  equal  to  R\  =  1  and  i?2  =  4,  respectively.  The  domain  is 
exactly  represented  by  a  single  NURBS  patch. 

The  whole  domain  boundary  is  assumed  to  be  clamped  and  we  assign  a  man¬ 
ufactured  solution  in  terms  of  displacement  components  (in  the  following  in¬ 
dicated  as  exact),  reading: 


u  =  (x2  +  y2  —  l)(x2  +  y2  —  16)  sin(a;)  sin (y), 
v  =  (x2  +  y2  —  l)(x2  +  y2  —  16)  sin(a;)  sin (y). 


(39) 


The  manufactured  solution  satisfies  the  prescribed  boundary  conditions  and 
the  load  is  computed  starting  from  the  manufactured  solution  and  imposing 
equilibrium. 


In  Figure  6  we  present  the  results  in  terms  of  relative  solution  error  in  the 
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Fig.  5.  Clamped  quarter  of  an  annulus.  Problem  geometry. 

L2-norm  versus  the  square  root  of  the  total  number  of  control  points  used, 
reporting  as  a  reference  the  convergence  rates  discussed  in  [3],  i.e. ,  p  and  p—  1 
for  even  and  odd  degree  p,  respectively. 


Fig.  6.  Clamped  quarter  of  an  annulus.  Error  plot  versus  the  square  root  of  number 
of  control  points  for  different  degree  NURBS. 


5.2  Single-material  single-patch  traction  test 


We  now  consider  a  square  domain  of  side  L  —  1,  subjected  to  uniform 
traction,  as  shown  in  Figure  7.  Accordingly,  we  specify  the  following  boundary 
conditions 

u  —  0  for  x  =  0  and  v  =  0  for  y  —  0, 
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while  we  assume  a  uniform  traction  q  for  x  =  L  and  traction-free  conditions 
for  y  =  L.  The  domain  consists  of  a  single  material  and  it  is  represented  by  a 
single  NURBS  patch. 

Such  a  problem  allows  us  to  test  the  proposed  numerical  scheme  for  Neumann 
boundary  conditions,  described  by  Equation  (27).  It  is  worth  emphasizing  that 
it  is  necessary  to  explicitly  impose  not  only  traction  boundary  conditions  (as 
in  classical  Galerkin  methods)  but  also  traction-free  ones  (which  are  instead 
naturally  satisfied  in  typical  Galerkin  methods).  Moreover,  the  problem  under 
investigation  is  also  characterized  by  a  corner  (point  A  in  Figure  7)  with  a 
combination  of  traction  boundary  condition  in  one  direction  and  traction-free 
in  the  other  direction,  which  corresponds  to  the  second  case  reported  in  Table 
1. 


Fig.  7.  Single-material  single-patch  traction  test.  Problem  geometry  and  boundary 
conditions. 


The  analytical  problem  solution  is  homogenous  and  governed  by  the  following 
set  of  equations 


\ 

r 

Gi 

1 
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0 

on 

f  22 

1 

1 

'v 

l-J 

1 

'v 

0 
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0~12 

(40) 


Being  in  a  plane-strain  situations,  enforcing  633  =  0,  it  is  possible  to  express 
033  in  terms  of  au  and  <r2 2  and  then,  requiring  <r22  =  oq2  =  0,  it  is  possible  to 
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compute  the  solution  as 


r  1  -  v2 
611  “ 

[  ^22  =  -  (1  +  Z')  (Til- 

Assuming  a  distributed  load  per  unit  length  q  =  10  and  material  constants 
E  =  1000  and  v  =  0.25,  the  displacement  components  of  point  A  are  then 

uA  =  9.375  x  1(T4  ,  vA  =  3.125  x  10“4. 

Such  an  analytical  solution  is  reproduced  up  to  machine  precision  by  the 
numerical  one  computed  using  a  single  element,  illustrating  the  good  behavior 
of  the  proposed  numerical  scheme  for  the  case  under  investigation. 

Figure  8  shows  the  horizontal  and  vertical  displacement  fields  (obtained  using 
p  =  q  =  2  and  3x3  control  points,  i.e. ,  one  element),  which  are  linear  in  the 
two  coordinate  variables,  as  expected. 


Fig.  8.  Single-material  single-patch  traction  test.  Horizontal  and  vertical  displace¬ 
ment  fields. 
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5.3  Two-material  two-patch  traction  test 


We  now  consider  a  rectangular  domain  fl,  as  sketched  in  Figure  9,  again 
subjected  to  uniform  traction.  However,  the  domain  is  now  assumed  to  consist 
of  two  material  subdomains. 
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Fig.  9.  Two-material  two-patch  traction  test.  Problem  geometry  and  boundary  con¬ 
ditions. 

The  problem  under  investigation  is  of  interest  with  respect  to  the  previous 
example,  since  it  introduces  a  boundary  point  (point  B  in  Figure  9)  with  a 
combination  of  a  traction-free  boundary  condition  and  an  interface  between 
different  materials. 

Similarly  to  what  has  been  done  in  the  previous  example,  the  idea  is  to  repro¬ 
duce  a  solution  homogeneous  in  the  y  direction  (and  piece-wise  homogeneous 
in  the  x  direction),  such  that  the  numerical  results  should  be  able  to  exactly 
reproduce  the  analytical  solution.  However,  to  obtain  such  a  solution,  it  is  nec¬ 
essary  to  properly  calibrate  the  elastic  constants.  Accordingly,  recalling  that 
Equations  40  are  valid  for  each  material,  we  require  the  transverse  strain  (i.e., 
the  strain  in  the  transverse  direction  with  respect  to  the  traction  direction)  to 
be  the  same  in  both  materials,  obtaining  the  following  relation 

Ei  _  ui  (1  +  zq) 

E2  v-2  (1  +  v 2)  ’ 

where  the  subscripts  indicate  the  material  numbers. 

For  the  problem  under  investigation  we  assume  U\  =  0.2,  v2  =  0.25  and  E2  = 
1000,  resulting  in  E\  =  768.  With  these  material  properties,  the  displacement 
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of  point  A  (indicated  in  Figure  9)  is  then 

uA  =  2.1875  x  1(T3,  vA  =  3.125  x  1(T4. 


We  solve  the  problem  numerically,  using  two  NURBS  patches  (i.e.,  a  patch 
for  each  material).  The  analytical  solution  is  matched  up  to  machine  precision 
by  the  numerical  one  computed  using  a  single  element  per  patch,  illustrat¬ 
ing  the  good  behavior  of  the  proposed  numerical  scheme  for  the  case  under 
investigation. 

Figure  10  shows  the  horizontal  and  vertical  displacement  fields  (obtained  using 
p  =  q  =  2  and  3x3  control  points  per  patch,  i.e.,  one  element  per  patch),  which 
are  linear  in  the  two  coordinate  variables  within  each  material,  as  expected.  We 
also  highlight  that,  as  desired,  a  perfectly  homogeneous  solution  is  obtained 
in  the  y  direction. 

We  also  notice  that  the  management  of  a  conforming  multi-patch  situation 
is  very  simple  in  the  proposed  collocation  method,  since  it  is  based  on  con¬ 
structing  the  discrete  equilibrium  relation  for  each  patch  and,  then,  summing 
the  equations  associated  to  collocation  points  shared  by  multiple  patches. 


5-4  Pressurized  thick-walled  cylinder  test 


We  now  consider  an  infinitely  long  and  internally  pressurized  thick-walled 
cylinder.  We  take  advantage  of  the  symmetry,  considering  only  a  quarter  of  the 
cylinder,  reducing  to  the  geometry  of  Figure  5  under  plane-strain  condition. 
We  impose  the  following  boundary  conditions 

u  —  0  for  x  —  0  and  v  —  0  for  y  —  0, 
and  assume  a  radial  pressure  load  P,  uniformly  distributed  at  the  inner  radius. 


For  the  problem  under  investigation  the  exact  solution  in  terms  of  radial  dis¬ 
placement  is 


PFP 


ur(r  = 


E  (Rl  -  Rj) 


(1  -  v)  r  T  (1  T  v)  — 
r 


where  r  is  the  radial  coordinate,  R{  and  R0  are  the  internal  and  the  outer 
radii,  and  E  and  v  are  the  Young’s  modulus  and  Poisson’s  ratio. 


Setting  Rt  =  1  and  R0  —  4,  E  —  1  and  v  =  0,  the  solution  becomes 
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Fig.  10.  Two-material  two-patch  traction  test.  Horizontal  and  vertical  displacement 
fields. 


Imposing  P  =  15/8,  results  in 


Mr  (1) 


17 


=  2.125 


ur( 4)  =  1. 


In  Figure  11  we  plot  the  displacement  magnitude  obtained  using  40  x  40 
control  points  and  p  =  q  =  4;  note  the  point-wise  (up  to  machine  precision) 
circumferential  symmetry  of  the  solution.  In  Figure  12  we  present  the  relative 
solution  error  in  the  L2-norm  versus  the  square  root  of  the  total  number  of 
control  points. 
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Fig.  11.  Pressurized  thick- walled  cylinder  test.  Displacement  magnitude. 


Fig.  12.  Pressurized  thick-walled  cylinder  test.  Error  plot  versus  the  square  root  of 
number  of  control  points  for  different  degree  NURBS. 


6  Numerical  tests:  dynamics 


To  solve  elastodynamic  problems,  we  employ  the  explicit  predictor  multi¬ 
corrector  algorithm  discussed  in  Section  4.  In  particular,  we  present  the  nu¬ 
merical  results  of  the  following  two  tests: 

•  a  clamped  rod,  excited  by  an  initial  velocity  distribution; 

•  a  clamped  plane-strain  quarter  of  an  annulus,  excited  by  a  time-dependent 
body  load  and  an  initial  velocity  distribution. 
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6. 1  Dynamics  of  a  clamped  rod 


We  consider  the  ID  elasto-dynamic  problem  of  a  rod  in  a  domain  [0, 1],  gov¬ 
erned  by  the  wave  equation 

il(x,t)  —  u"(x,t)  —  0  \/x  e  (0, 1),  Vt  e  (0, T),  (41) 

where  u  =  u(x,t)  is  the  unknown  displacement  in  terms  of  the  axial  coor¬ 
dinate  x  and  time  t,  while  (•)  and  (•)'  represent  time  and  space  derivatives, 
respectively. 

We  consider  the  following  initial  conditions 

u(x,  0)  =  0,  ii(x,  0)  =  27r  sin(27nr),  Vo;  G  (0,1),  (42) 

and  boundary  conditions 

u(0,t)  =  u{l,t)  =  0,  Vf  G  (0,T),  (43) 

such  that  the  exact  solution  is 

u(x,  t )  =  sin(27rx)  sin(27rf).  (44) 


In  Figure  13  we  present  the  relative  solution  error  in  the  L2-norm  at  the 
final  time  T  (T  =  7/4  in  our  simulations),  plotted  versus  the  total  number 
of  control  points.  Such  numerical  solutions  are  obtained  using  the  explicit 
predictor  multi-corrector  algorithm  with  two  corrector  passes.  We  remark  that, 
since  the  adopted  algorithm  is  only  second-order  in  time,  in  order  to  manifest 
the  expected  high-order  convergence  rates  in  space,  the  time  step  and  spatial 
mesh-size  must  be  selected  such  that  the  space  discretization  errors  dominate 
the  time  discretization  errors. 


6.2  Dynamics  of  a  clamped  plane- strain  quarter  of  an  annulus 


We  consider  a  2D  plane-strain  elastodynamic  problem  where  a  clamped  quar¬ 
ter  of  an  annulus  (whose  geometry  is  described  in  Figure  5)  is  excited  by  a 
time-dependent  body  force  and  an  initial  velocity  distribution  such  that  the 
exact  solution  is 


u  =  (x2  +  y2  —  l)(x2  +  y2  —  16)  sin(x)  sin (y)  sin(27r t), 
v  =  (x2  +  y2  —  l)(x2  +  y2  —  16)  sin(x)  sin (y)  sin(27rf). 


(45) 
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Fig.  13.  Dynamics  of  a  clamped  rod.  Error  plot  versus  number  of  control  points  for 
different  degree  NURBS. 


In  Figure  14  we  present  the  relative  solution  error  in  the  L2-norm  at  the  final 
time  T  (T  =  7/4  in  our  simulations),  plotted  versus  the  square  root  of  the  total 
number  of  control  points.  Again,  such  numerical  solutions  are  obtained  using 
the  explicit  predictor  multi-corrector  algorithm  with  two  corrector  passes.  The 
numerical  results  are  consistent  with  those  obtained  in  the  ID  case.  Also 
here,  the  time  step  and  spatial  mesh-size  are  selected  in  order  that  the  space 
discretization  errors  dominate  the  time  discretization  errors. 


Fig.  14.  Dynamics  of  a  clamped  plane-strain  quarter  of  an  annulus.  Error  plot  versus 
the  square  root  of  number  of  control  points  for  different  degree  NURBS. 
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7  Conclusions 


In  this  paper  we  have  developed  IGA  collocation  formulations  applicable 
to  multi-patch  NURBS  configurations,  investigated  the  treatment  of  various 
boundary  and  patch  interface  conditions,  and  extended  the  methodology  to 
explicit  dynamics.  We  have  argued  that  the  procedures  are  stable,  robust, 
higher-order  accurate  and  efficient  in  the  sense  that  they  involve  a  minimum 
number  of  quadrature  points,  and  we  have  presented  numerical  results  to  sup¬ 
port  our  claims. 

We  believe  that  the  new  methods  have  the  potential  to  offer  a  superior  alterna¬ 
tive  to  existing  finite  element  technology  based  on  one-point  Gaussian  quadra¬ 
ture  bilinear  quadrilateral  and  trilinear  hexahedral  elements.  However,  much 
research  still  remains  to  be  done.  We  believe  that  the  focus  of  this  work  should 
be  on  fully  nonlinear  problems,  shell  formulations,  three-dimensional  solids, 
industrial  scale  calculations,  and  extension  to  hierarchically  refined  NURBS 
and  T-splines,  which  we  hope  to  pursue  in  the  future. 


A  Appendix:  Construction  of  Sfjk 


We  do  not  show  the  proof  of  the  following  classical  result. 


Lemma  A.l  Let  f,g  be  real  valued  functions  defined  on  a  neighborhood  of 
Tij'k  in  fb  Let  f  be  integrable  and  continuous  in  T%.hk,  and  let  g  be  in  L°°. 


Then 


lim 

€ — >0 


^supp(r) 

[  g  STk  <C  V £  6  (0, 1]  , 

./supple;) 


(A.l) 


with  C  a  constant  independent  of  e. 


Moreover,  the  following  result  is  easy  to  prove. 

Lemma  A. 2  Let  7  be  a  C1  curve  originating  at  Tl^  k.  Let  f  be  an  integrable 
real  valued  function  defined  on  the  curve,  continuous  in  Then 

£  So  /,  f  5 ’*j’k  =  Cs  f(Ti’Lk)  ’  (A-2) 

with  C5  a  constant  independent  of  /,  £,7. 
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